#  Clear memory
rm(list=ls())

# Load packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  tidyverse, data.table)

# Import fips code list 
fips_order = fread("raw/modified-ap3/fips.csv") 
fips_order = fips_order %>% 
  mutate(order = seq(1:nrow(fips_order))) %>%
  rename(fips = V1)

# Import counterfactual to get fips order for subset of counties we need
counterfactual = fread("raw/modified-ap3/counterfactual_fips.csv")

# Get MD from inverted AP3 run for 1997.

nox_md = fread("data/modified-ap3/NOX_md_L_1997_prime_aged.csv")  %>% 
  mutate(order = seq(1:nrow(fips_order))) %>%
  left_join(fips_order, by = c("order")) %>%
  select(fips, starts_with("V")) %>%
  rename(fips_origin = fips) %>%
  pivot_longer(cols = V1:V3109, names_to = "name", values_to = "damage") %>%
  mutate(order = rep(seq(1:nrow(fips_order)), times = 3109)) %>%
  left_join(fips_order, by = c("order")) %>%
  rename(fips_destination = fips) %>%
  select(fips_origin, fips_destination, damage) %>%
  rename(fips = fips_destination) %>%
  right_join(counterfactual) %>%
  pivot_wider(names_from = fips, values_from = damage, names_prefix = "V_") %>%
  rename(fips = fips_origin) %>%
  right_join(counterfactual) %>%
  select(-c(fips))

fwrite(nox_md, file = "data/modified-ap3/NOX_md_L_1997_prime_aged_for_counterfactual.csv")

nh3_md = fread("data/modified-ap3/NH3_md_L_1997_prime_aged.csv") %>% 
  mutate(order = seq(1:nrow(fips_order))) %>%
  left_join(fips_order, by = c("order")) %>%
  select(fips, starts_with("V")) %>%
  rename(fips_origin = fips) %>%
  pivot_longer(cols = V1:V3109, names_to = "name", values_to = "damage") %>%
  mutate(order = rep(seq(1:nrow(fips_order)), times = 3109)) %>%
  left_join(fips_order, by = c("order")) %>%
  rename(fips_destination = fips) %>%
  select(fips_origin, fips_destination, damage) %>%
  rename(fips = fips_destination) %>%
  right_join(counterfactual) %>%
  pivot_wider(names_from = fips, values_from = damage, names_prefix = "V_") %>%
  rename(fips = fips_origin) %>%
  right_join(counterfactual) %>%
  select(-c(fips))

fwrite(nh3_md, file = "data/modified-ap3/NH3_md_L_1997_prime_aged_for_counterfactual.csv")

pm25_md = fread("data/modified-ap3/PM25_md_L_1997_prime_aged.csv") %>% 
  mutate(order = seq(1:nrow(fips_order))) %>%
  left_join(fips_order, by = c("order")) %>%
  select(fips, starts_with("V")) %>%
  rename(fips_origin = fips) %>%
  pivot_longer(cols = V1:V3109, names_to = "name", values_to = "damage") %>%
  mutate(order = rep(seq(1:nrow(fips_order)), times = 3109)) %>%
  left_join(fips_order, by = c("order")) %>%
  rename(fips_destination = fips) %>%
  select(fips_origin, fips_destination, damage) %>%
  rename(fips = fips_destination) %>%
  right_join(counterfactual) %>%
  pivot_wider(names_from = fips, values_from = damage, names_prefix = "V_") %>%
  rename(fips = fips_origin) %>%
  right_join(counterfactual) %>%
  select(-c(fips))

fwrite(pm25_md, file = "data/modified-ap3/PM25_md_L_1997_prime_aged_for_counterfactual.csv")

so2_md = fread("data/modified-ap3/SO2_md_L_1997_prime_aged.csv") %>% 
  mutate(order = seq(1:nrow(fips_order))) %>%
  left_join(fips_order, by = c("order")) %>%
  select(fips, starts_with("V")) %>%
  rename(fips_origin = fips) %>%
  pivot_longer(cols = V1:V3109, names_to = "name", values_to = "damage") %>%
  mutate(order = rep(seq(1:nrow(fips_order)), times = 3109)) %>%
  left_join(fips_order, by = c("order")) %>%
  rename(fips_destination = fips) %>%
  select(fips_origin, fips_destination, damage) %>%
  rename(fips = fips_destination) %>%
  right_join(counterfactual) %>%
  pivot_wider(names_from = fips, values_from = damage, names_prefix = "V_") %>%
  rename(fips = fips_origin) %>%
  right_join(counterfactual) %>%
  select(-c(fips))

fwrite(so2_md, file = "data/modified-ap3/SO2_md_L_1997_prime_aged_for_counterfactual.csv")

voc_md = fread("data/modified-ap3/VOC_md_L_1997_prime_aged.csv") %>% 
  mutate(order = seq(1:nrow(fips_order))) %>%
  left_join(fips_order, by = c("order")) %>%
  select(fips, starts_with("V")) %>%
  rename(fips_origin = fips) %>%
  pivot_longer(cols = V1:V3109, names_to = "name", values_to = "damage") %>%
  mutate(order = rep(seq(1:nrow(fips_order)), times = 3109)) %>%
  left_join(fips_order, by = c("order")) %>%
  rename(fips_destination = fips) %>%
  select(fips_origin, fips_destination, damage) %>%
  rename(fips = fips_destination) %>%
  right_join(counterfactual) %>%
  pivot_wider(names_from = fips, values_from = damage, names_prefix = "V_") %>%
  rename(fips = fips_origin) %>%
  right_join(counterfactual) %>%
  select(-c(fips))

fwrite(voc_md, file = "data/modified-ap3/VOC_md_L_1997_prime_aged_for_counterfactual.csv")

nox_transport = fread("data/modified-ap3/NOX_transport_L_1997_prime_aged.csv")  %>% 
  mutate(order = seq(1:nrow(fips_order))) %>%
  left_join(fips_order, by = c("order")) %>%
  select(fips, starts_with("V")) %>%
  rename(fips_origin = fips) %>%
  pivot_longer(cols = V1:V3109, names_to = "name", values_to = "damage") %>%
  mutate(order = rep(seq(1:nrow(fips_order)), times = 3109)) %>%
  left_join(fips_order, by = c("order")) %>%
  rename(fips_destination = fips) %>%
  select(fips_origin, fips_destination, damage) %>%
  rename(fips = fips_destination) %>%
  right_join(counterfactual) %>%
  pivot_wider(names_from = fips, values_from = damage, names_prefix = "V_") %>%
  rename(fips = fips_origin) %>%
  right_join(counterfactual) %>%
  select(-c(fips))

fwrite(nox_transport, file = "data/modified-ap3/NOX_transport_L_1997_prime_aged_for_counterfactual.csv")

nh3_transport = fread("data/modified-ap3/NH3_transport_L_1997_prime_aged.csv") %>% 
  mutate(order = seq(1:nrow(fips_order))) %>%
  left_join(fips_order, by = c("order")) %>%
  select(fips, starts_with("V")) %>%
  rename(fips_origin = fips) %>%
  pivot_longer(cols = V1:V3109, names_to = "name", values_to = "damage") %>%
  mutate(order = rep(seq(1:nrow(fips_order)), times = 3109)) %>%
  left_join(fips_order, by = c("order")) %>%
  rename(fips_destination = fips) %>%
  select(fips_origin, fips_destination, damage) %>%
  rename(fips = fips_destination) %>%
  right_join(counterfactual) %>%
  pivot_wider(names_from = fips, values_from = damage, names_prefix = "V_") %>%
  rename(fips = fips_origin) %>%
  right_join(counterfactual) %>%
  select(-c(fips))

fwrite(nh3_transport, file = "data/modified-ap3/NH3_transport_L_1997_prime_aged_for_counterfactual.csv")

pm25_transport = fread("data/modified-ap3/PM25_transport_L_1997_prime_aged.csv") %>% 
  mutate(order = seq(1:nrow(fips_order))) %>%
  left_join(fips_order, by = c("order")) %>%
  select(fips, starts_with("V")) %>%
  rename(fips_origin = fips) %>%
  pivot_longer(cols = V1:V3109, names_to = "name", values_to = "damage") %>%
  mutate(order = rep(seq(1:nrow(fips_order)), times = 3109)) %>%
  left_join(fips_order, by = c("order")) %>%
  rename(fips_destination = fips) %>%
  select(fips_origin, fips_destination, damage) %>%
  rename(fips = fips_destination) %>%
  right_join(counterfactual) %>%
  pivot_wider(names_from = fips, values_from = damage, names_prefix = "V_") %>%
  rename(fips = fips_origin) %>%
  right_join(counterfactual) %>%
  select(-c(fips))

fwrite(pm25_transport, file = "data/modified-ap3/PM25_transport_L_1997_prime_aged_for_counterfactual.csv")

so2_transport = fread("data/modified-ap3/SO2_transport_L_1997_prime_aged.csv") %>% 
  mutate(order = seq(1:nrow(fips_order))) %>%
  left_join(fips_order, by = c("order")) %>%
  select(fips, starts_with("V")) %>%
  rename(fips_origin = fips) %>%
  pivot_longer(cols = V1:V3109, names_to = "name", values_to = "damage") %>%
  mutate(order = rep(seq(1:nrow(fips_order)), times = 3109)) %>%
  left_join(fips_order, by = c("order")) %>%
  rename(fips_destination = fips) %>%
  select(fips_origin, fips_destination, damage) %>%
  rename(fips = fips_destination) %>%
  right_join(counterfactual) %>%
  pivot_wider(names_from = fips, values_from = damage, names_prefix = "V_") %>%
  rename(fips = fips_origin) %>%
  right_join(counterfactual) %>%
  select(-c(fips))

fwrite(so2_transport, file = "data/modified-ap3/SO2_transport_L_1997_prime_aged_for_counterfactual.csv")

voc_transport = fread("data/modified-ap3/VOC_transport_L_1997_prime_aged.csv") %>% 
  mutate(order = seq(1:nrow(fips_order))) %>%
  left_join(fips_order, by = c("order")) %>%
  select(fips, starts_with("V")) %>%
  rename(fips_origin = fips) %>%
  pivot_longer(cols = V1:V3109, names_to = "name", values_to = "damage") %>%
  mutate(order = rep(seq(1:nrow(fips_order)), times = 3109)) %>%
  left_join(fips_order, by = c("order")) %>%
  rename(fips_destination = fips) %>%
  select(fips_origin, fips_destination, damage) %>%
  rename(fips = fips_destination) %>%
  right_join(counterfactual) %>%
  pivot_wider(names_from = fips, values_from = damage, names_prefix = "V_") %>%
  rename(fips = fips_origin) %>%
  right_join(counterfactual) %>%
  select(-c(fips))

fwrite(voc_transport, file = "data/modified-ap3/VOC_transport_L_1997_prime_aged_for_counterfactual.csv")
